Identifying Genes Associated with Female Flower Development of Phellodendron amurense Rupr. Using a Transcriptomics Approach

Phellodendron amurense Rupr., a species of Rutaceae, is a nationally protected and valuable medicinal plant. It is generally considered to be dioecious. With the discovery of monoecious P. amurense, the phenomenon that its sex development is regulated by epigenetics has been revealed, but the way epigenetics affects the sex differentiation of P. amurense is still unclear. In this study, we investigated the effect of DNA methylation on the sexual development of P. amurense. The young inflorescences of male plants were treated with the demethylation agent 5-azaC, and the induced female flowers were obtained. The induced female flowers’ morphological functions and transcriptome levels were close to those of normally developed plants. Genes associated with the development of female flowers were studied by comparing the differences in transcriptome levels between the male and female flowers. Referring to sex-related genes reported in other plants, 188 candidate genes related to the development of female flowers were obtained, including sex-regulating genes, genes related to the formation and development of sexual organs, genes related to biochemical pathways, and hormone-related genes. RPP0W, PAL3, MCM2, MCM6, SUP, PIN1, AINTEGUMENTA, AINTEGUMENTA-LIKE6, AGL11, SEUSS, SHI-RELATED SEQUENCE 5, and ESR2 were preliminarily considered the key genes for female flower development. This study has demonstrated that epigenetics was involved in the sex regulation of P. amurense, with DNA methylation as one of its regulatory modes. Moreover, some candidate genes related to the sexual differentiation of P. amurense were obtained with analysis. These results are of great significance for further exploring the mechanism of sex differentiation of P. amurense and studying of sex differentiation of plants.


Introduction
As the sex organ of plants, floral organs play an important role in plant reproduction. Flower development begins with flower bud differentiation. During the development of flower buds, the meristems differentiate into different flower organ primordia first and then develop into corresponding flower organs. Flower buds can be divided into bisexual flowers and unisexual flowers according to the development of pistils and stamens, and unisexual flowers can be divided into female and male flowers according to the abortion of stamens or pistils. Genetic factors, environmental conditions, hormones, and other factors can regulate the sex differentiation in flowers. Among them, epigenetics has been widely considered to play a crucial role in plant sexual differentiation. DNA methylation is one of the most intensively studied modes of epigenetic regulation [1,2].
DNA methylation is a universal regulatory mode in plants and animals. DNA methylation is the process of transferring methyl groups to a DNA nucleotide catalyzed by DNA methyltransferase with S-adenosylmethionine as methyl donor, which can always influence the expression of target genes [3]. DNA methylation can regulate plant sex by modifying sex-determining genes. As the sex switch in poplar, the ARR17 gene is silenced by malespecific DNA methylation [4]; the transition from male to female flowers in gynoecious melon results from DNA methylation in the promoter of CmWIP1 that can contribute to the development of male flowers through carpel abortion [5]; and the male flowers of Diospyros kaki are closely connected with the low expression of the female-determining gene MeGI because of DNA methylation [6,7]. In addition, DNA methylation can regulate plant sex by modifying sex-related genes in response to external cues. DNA methylation of CpHUA1 in papaya in response to temperature-related stress can lead to sex reversal from male to hermaphrodite [8]. The sex of cucumber is unstable and easily influenced by temperature and photoperiod. DNA methylation plays a key role in regulating cucumber sex by affecting sex-related genes in response to external conditions [9,10]. With the development of research technology, an increasing number of plants have been reported to regulate their sex with epigenetics.
P. amurense Rupr., a lofty deciduous arbor in the Rutaceae, is a valuable medicinal tree species with high economic value [11]. It is recorded that P. amurense is a dioecious plant with terminal inflorescences. The stamens or pistils degenerate during the development of unisexual flowers [12]. However, monoecious P. amurense was first discovered by Fan et al. [13]. The top of monoecious P. amurense died when it was young, and the lateral buds grew and developed into branches of different sexes, which suggested that the sex formation of P. amurense occurred after the formation of branches, and there was different sex expression in the same genetic background. The sex expression of P. amurense is closely related to epigenetic regulation. DNA methylation is a common epigenetic regulatory mode that has been extensively studied. Therefore, we aimed to explore the sexual differentiation in P. amurense through demethylation in this study.
The common demethylation reagent 5-azacytidine (5-azaC) can specifically inhibit DNA methyltransferase from preventing DNA methylation [14,15]. Many studies have shown that 5-azaC can significantly reduce the methylation level in plants and affect their growth and development [16][17][18]. This method is also gradually applied in the study of plant sexual differentiation. Research on Melandrium album has shown that treatment with 5-azaC could induce a sex change to andromonoecy in about 21% of male plants, while no apparent phenotypic effect was observed in females [19]. In this study, we treated male flowers of P. amurense with 5-azaC and found that some of these male flowers had undergone sex transition, forming 5-azaC-induced female flowers.
In this study, the generation of 5-azaC-induced female flowers suggests that the sex of P. amurense is regulated by epigenetics. To explore genes related to the development of female flowers, we have compared the transcriptome data of female flowers and induced female flowers with male flowers, in combination with the sex-related genes reported in the literature. The results of this study can lay a foundation for further exploring the sexual differentiation mechanism of P. amurense and provide important information for the study of plant sexual differentiation.

Plant Materials and 5-azaC Treatment
The flower buds of 11 male P. amurense individuals preserved in Beijing Medical Botanical Garden were treated with 1 mM 5-azaC in sterile distilled water (floral spray of about 0.3 mL 5-azaC to every inflorescence), which continued during the dormancy period (16 November 2021-10 January 2022) and the growth period (8 March 2022-20 April 2022) once every 2 days. Twenty-eight and twenty-two treatments were performed in the dormant period and the growing period, respectively.

Samples Collection
When the sex of flowers could be accurately distinguished in appearance (diameter of flower buds of about 5 mm), induced female flowers (F i ) and male flowers without sex conversion after treatment (M 5-azaC ) were collected, and male flowers without 5-azaC treatment on the same plant were collected as controls (male flowers for CK; M ck ). Female flowers (F) and male flowers (M) of nontreated plants were collected simultaneously. Each sample had three biological repeats. The samples above were immediately frozen in liquid nitrogen and stored at −80 • C.

RNA Isolation, Illumina Sequencing, Transcriptome Assembly, and Annotation
Total RNAs were extracted with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. The quality and quantity of RNA samples were assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).
Total RNA was used as input material for the RNA sample preparations. Poly (A) mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Then, a cDNA library was constructed. First-strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase, then using RNaseH to degrade the RNA. Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and dNTPs. After passing the library quality inspection, the library preparations were sequenced on an Illumina Novaseq 6000 platform (Illumina, San Diego, CA, USA), and 150 bp paired-end reads were generated. Clean data (clean reads) with high quality were obtained by removing reads containing adapters, reads containing N base, and low-quality reads from raw data. Transcriptome assembly was accomplished using Trinity v2.4.0 [20]. The longest non-redundant unigenes were acquired by removing sequence splicing and redundancy using Corset v4.6 [21]. The assembly quality was assessed with BUSCO [22].
For further annotation of unigenes, all the assembly unigenes were first searched against the NCBI non-redundant protein sequences (Nr) using diamond v0.8.22 (e-value < e −5 ) [23]. Then, Blast2GO v2.5 [24] was used to obtain the Gene Ontology (GO) annotation of unigenes based on the Nr annotation (e-value < e −6 ). Pathway assignments were carried out according to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database using the KEGG Automatic Annotation Server (e-value < e −10 ) [25].

Acquisition of Sex-Related Genes in P. amurense
It has been shown that the underlying mechanisms controlling flower development are largely conserved in distantly related dicotyledonous plant species [26]. Therefore, genomic resources generated from other plants could be used to identify the potential genes involved in the sexual differentiation and flower development of P. amurense. A literature survey was undertaken to list 1281 genes involved in sex differentiation and flower development in other plants, involving 189 papers from more than 30 journals (the latest published in 2022) [4,8,. These can be divided into the following four categories according to function: sex regulatory genes (Group A), including sex-determining genes, genes involved in sex regulatory mechanisms, and genes located in sex-determining regions (SDR); genes related to floral initiation and development (Group B), such as genes related to floral primordium formation and floral organ development; genes related to biochemical metabolic pathways (Group C); and genes related to plant hormones (Group D) (Table S1).
Nucleotide sequences of the above sex-related genes reported in the literature were downloaded from the NCBI Genbank database in FASTA format. A series of BLASTN analyses using e-value < e −10 as a threshold identified broadly conserved sequences of potential sex-related genes from the P. amurense transcriptome [216]. The sequences acquired with BLASTN were input to the Conserved Domains platform of NCBI to ensure the homologous sequences shared the same domains.

Differential Expression Analysis
RSEM was used for transcript abundance estimation [217]. Pre-processed RNA-Seq paired reads for each flower sex type were mapped to the final assembled transcriptome using Trinity. All read counts were calculated using the fragments per kilobase of transcript per million fragments mapped (FPKM) method. Differential expression analysis was conducted between female and male flowers in nontreated trees. To obtain information on sex-related genes of the female and male P. amurense, differentially expressed genes (DEGs) between female and male flowers were identified by comparing F with M using DESeq2 [218] based on criteria set as a |log2 fold change| ≥ 1 and a false discovery rate (FDR) ≤ 0.05. Further, we mainly aimed at the sex-related genes identified in P. amurense transcriptome using blasting to screen genes in a small scope. Therefore, we conducted a differential expression analysis of sex-related genes between F and M using FDR ≤ 0.05 as the criterion.
An analysis of genes related to female flower development was carried out. First, differential expression analysis of sex-related genes between F i and M ck was conducted to obtain sex-related genes that changed during sex conversion after treating male flowers with 5-azaC. The FDR ≤ 0.05 was used as a threshold. Next, considering that some genes unrelated to sex might also be changed with 5-azaC, those genes needed to be removed. Because there was no sex change in part of male flowers treated with 5-azaC (M 5-azaC ), DEGs between M 5-azaC and M ck could be regarded to be unrelated to sexual differentiation. After removing DEGs between M 5-azaC and M ck from DEGs between F i and M ck , the remaining genes were directly related to female flower development. The key genes related to female flower development were explored in combination with gene functions.

Quantitative Real-Time PCR (qRT-PCR) Analysis
Thirty-six genes were selected to validate their expression patterns. qPCR primers were designed by prime primer 5 (Table S2) [219]. The total RNA of each sample was reverse-transcripted into cDNA using PrimeScript™ RT reagent Kit with gDNA Eraser (Takara, Dalian, China). Real-time qPCR was performed using TB Green ® Premix Ex Taq™ (Tli RNaseH Plus; Takara, Dalian, China) and analyzed on a BIORAD CFX Real-Time System. Ubiquitin was used for normalization, and the expression ratio was calculated using the 2 −∆∆Ct formula.

Generation of Induced Female Flowers
After treatment with 5-azaC, some male flower buds were induced to transform into female flowers with complete structure and function, while the remaining male flowers had no change in structure and function, and there were no other variation types. Induced female flowers were generated in about 70% of male individuals treated, and induced female flower buds were generated in 7-50% of inflorescence treated among the male individuals in which induced female flowers were generated. Induced female flowers were completely identified as natural female flowers in morphology, structure, and function, characterized in that stamens degenerated to different degrees and ovaries could develop into fruits with seeds ( Figure 1).

Results of Transcriptome Assembly and Annotation
The 15 cDNA libraries constructed from all samples were sequenced on an Illumina high-throughput sequencing platform. After filtering the adaptors and low-quality sequences, 335,043,541 high-quality clean reads were generated from the 15 cDNA libraries. The resulting assembly consisted of 263,986 transcripts with an N50 value of 2180 bp. Removing redundant sequences resulted in a total of 86,610 unigenes with an N50 value of 1790 bp. The length distribution of unigenes ranged from 301 to 16,091 bp, nearly 33% of which were 300-500 bp long, approximately 50% of which were also distributed between 500 and 2000 bp, while the remaining 17% were clustered under a size distribution of ≥2000 bp.
A total of 52,231 unigenes (60.30%) were annotated in at least one of the Nr, GO, and KEGG databases ( Figure 2). There were 47,588 (54.94%), 31,006 (35.79%), and 17,184 (19.84%) unigenes annotated in the Nr, GO, and KEGG databases, respectively. According to the species distribution analysis of the Nr database, the top hit was from Citrus sinensis (25.4%). Unigenes were also matched significantly with Citrus clementina (20.8%) and Citrus unshiu (20.4%); additionally, a small number of unigenes matched with Vitis vinifera (3.0%) and Quercus variabilis (1.6%). The result of Nr annotation suggested P. amurense was more similar to C. sinensis, Citrus clementina, and Citrus unshiu, which was consistent with the fact that those four species all belong to Rutaceae ( Figure 2).
The GO analysis categorized 31,006 unigenes into three main categories (biological processes, BP; cellular components, CC; and molecular functions, MF). Among the biological processes category, a total of 22,827 unigenes were categorized into 26 functional groups, and "cellular processes" (GO: 0009987) and "metabolic processes" (GO: 0008152) were the predominant groups. Within the cellular components category, a total of 25,399 unigenes were categorized into five functional groups, and cellular anatomical entry (GO: 0110165), intracellular (GO: 0005622), and protein-containing complex (GO: 0032991) were the most overrepresented groups. For the molecular function category, a total of 14,535 unigenes were categorized into 12 functional groups, and a large proportion of unigenes were clustered into binding (GO: 0005488) and catalytic activity (GO: 0003824) ( Figure 2).

Differential Expression Analysis of Transcriptome between Female and Male Flowers
Differential expression analysis was conducted between female and male flowers (F vs. M) to acquire information on DEGs between the two samples. A total of 16,600 DEGs (19.16%) were identified between F and M, including 7107 up-regulated genes and 9493 down-regulated genes in F compared with M.
The KEGG pathways in which DEGs enriched significantly (corrected p-value < 0.05) and the top 30 GO terms distribution for DEGs are shown in Figure 3. Red and green dots present female highly expressed and male highly expressed genes, respectively, and blue dots mean no significant difference. (B,C) Enrichment of female highly expressed and male highly expressed genes in KEGG pathways. (D,E) Enrichment of female highly expressed and male highly expressed genes in the GO database. * significant enrichment (corrected p-value < 0.05).

Identification of Sex-Related Genes in Female and Male Flowers and Differential Expression Analysis
To identify genes potentially involved in flower development and sex differentiation, the homologous sequences of sex-related genes from P. amurense were acquired using BLASTN (Table S3). A total of 604 and 584 homologous sequences were acquired from F and M, respectively. F and M shared similar distributions of several genes in the four gene types, with the largest proportion of Group B and the smallest proportion of Group A (Table 1). To further explore the genes related to flower development and sex differentiation in P. amurense, differential expression analysis of sex-related genes between F and M was carried out. Using an FDR of 0.05, a total of 332 sex-related DEGs were identified between F and M, including 172 up-regulated genes and 160 down-regulated genes in F compared with M. The log2FC values of the DEGs ranged from 9.23 to −19.12 (Table 2 and Table S4). Some DEGs with female or male-specific expressions might strongly contribute to sexual differentiation. The female-specific genes included early development genes of female organs ESR2 and SUP, chemical metabolism genes MYBC and ATB51, indoleacetic acid (IAA)-related genes IAA27 and SAU50, ethylene (ETH)-related gene AIL5, and gibberellin (GA)-related gene G3OX3. The male-specific genes included stamen development gene GNL2, pollen maturation and pollen tube growth gene AGL104, pollen early development gene ZAT3, male gametophyte development gene MMD1, pectinesterase-related gene PEI, and ATPase-related gene V-type proton ATPase subunit G1.
Among A-group DEGs, the sex-determining gene TOZ19 and genes involved in plant sex regulation, such as RPP0W, PAL3, MCM6, MCM2, and ABA2, had higher expression in F, while male-determining gene GPAT3 and the sex regulatory gene GAST1 had higher expression in M.

Identification of Key Genes Related to Female Flower Development
Some male flowers could be induced into female flowers with no other variation types after treatment with 5-azaC. Key genes related to female flower development in P. amurense might be involved in the sex conversion of male flowers. A total of 612, 611, and 606 homologous sequences of sex-related genes were acquired from F i , M ck, and M 5-azac , respectively, using BLASTN. The three samples shared similar distributions in the number of genes in the four gene types (Table 1).
Using an FDR of 0.05, 296 sex-related DEGs were identified between F i and M ck , including 142 up-regulated genes and 154 down-related genes in F i compared with M ck . The log2FC values of the DEGs ranged from 10.36 to −13.55 (Table 2 and Table S5).
F i was close to F at the transcriptome level; they shared 216 DEGs, which were regarded as candidate genes for female flower development compared with male flowers. They were more reliable genes for female flower development.
As a demethylation reagent, 5-azaC could provoke extensive physiological changes not limited to sex. Correlation analysis using the Spearman method and principal component analysis showed a good correlation between the replicate sets of F i and M ck , but there was an obvious difference among the three biological repeats of M 5-azaC , which were close to F i or M ck , respectively. That suggested that changes in varying degrees had happened within M 5-azaC after the treatment with 5-azaC, but these changes might involve other biological activities unrelated to sex formation ( Figure 4). Therefore, the DEGs between M 5-azaC and M ck were less likely to be key genes and could be regarded as not linked with sex formation. A total of 28 such genes were removed from 216 candidate genes of female flowers, and finally, 188 candidate genes were screened out ( Figure 5 and Table S6). The candidate genes were mainly involved in sex regulation, the initiation and development of the floral organ, biochemical metabolic pathways, and plant hormones. Combining with their functions, genes associated with sex regulation and floral initiation, including RPP0W, PAL3, MCM2, MCM6, SUP, ANT, AIL6, AGL11, SEUSS, SRS5, and ESR2, were preliminarily considered to be the key genes for female flower development.

Verification of DEGs Using qRT-PCR Analysis
We paid attention to the following two types of genes when selecting genes to be validated. One was A-group genes because the genes involved in sex regulation might contribute strongly to sex determination. Another was up-regulated genes in F i , considering that demethylation is more likely to improve the expression levels of target genes. Therefore, 36 genes, including A-group genes and up-regulated genes in F i , were selected for qRT-PCR validation ( Figure 6).

Discussion
RNA-Seq has been a common technique applied in plants for exploring DEGs between female and male flowers, such as papaya [220] and Cannabis sativa [68]. However, there have been few reports on genes related to female flower development in which demethylation was used to generate induced female flowers and further analysis was performed combining sex-related genes known in other plants. In this study, we succeeded at inducing male flowers into female flowers using the demethylation reagent 5-azaC. After acquiring sex-converted materials, 188 candidate genes of female flower development were identified with further analysis.
In P. amurense, the female flowers are present as bisexual tissue at the initial stage, and then the stamen stops developing, forming a unisexual flower. This process requires many specific genes to participate in each development stage. The 188 candidate genes acquired in this study are involved in sex regulation, flower development, biochemical pathway, and hormone regulation. The up-regulated genes are more likely to be influenced by demethylation. Among them, PAL3, RPP0W, MCM2, MCM6, and SRS5 might be involved in sexual differentiation regulation in P. amurense during the sex-determining stage. PAL3, encoding a key player in the phenylpropanoid pathway, was detected as a candidate gene for female flower development in Populus tomentosa with a sex-specific methylation alteration [221]. RPP0W, identified in the SDR (sex-determining region) cassette occurring in nearly all females and never in males, potentially regulated the sexual differentiation in Fragaria [162]. MCM2 and MCM6, two genes involved in the cell cycle pathway, had much higher expression in the gynoecious line than those in weak female cucumbers and were considered key candidate genes of sexual differentiation in cucumbers [166]. SRS5 had higher expression in early gynoecious inflorescence buds than monoecious plants in Jatropha curcas, which was considered a candidate regulator of the sex determination [222].
During the flowering stage in P. amurense, ANT, AIL6, SEUSS, AGL11, SUP, PIN1, and ESR2 might promote the initiation of female floral organ primordia. ANT and its paralog AIL6 are two key regulators of flower development, including floral organ initiation, identity specification, growth, and patterning. On the one hand, they are important for establishing the flower primordia, such as ovule and female gametophyte in Arabidopsis. On the other hand, they participate in subsequent female flower development by regulating target genes including SPATULA, YABBY, SEP3, PHB, AS1, REV, and FIL and genes associated with hormones identified in P. amurense [223][224][225]. SEUSS cooperates with ANT in a partially redundant manner to regulate the expression of downstream genes critical for the formation of ovules [226]. AGL11 participates in the early development of ovules, which was demonstrated in many species, including Arabidopsis thaliana and Punica granatum [227,228]. SUP and PIN1 were demonstrated to be involved with reproductive phase transition and female flower transition by promoting the abortion of male flower primordia in Jatropha curcas [229]. Furthermore, SUP plays a role in maintaining the boundaries between stamens and carpels and regulating the development of outer ovule integument in Arabidopsis [230]. ESR2 is a member of transcription factor BOL/DRNL/ESR2/SOB, which is expressed at the very early stages in aerial organ formation and has been proposed to be a marker for organ founder cells [231]. The results of this study suggest that the genes above were likely to participate in the early establishment of female flowers, which have been regarded as key genes for female flower development and likely make significant contributions to female formation in P. amurense. It is worth mentioning that PAL3, MCM2, MCM6, AGL11, ANT, AIL6, SRS5, ESR2, and SUP have shown differential expression during our research aimed at early flowers (the result is unpublished), which could make the results of this study more convincing to some degree. Moreover, CYC2CL, RL1, RL2, PROLIFERA, FBP2, YABBY1, YABBY2, FIL, AS1, and PHB might be involved in subsequent female flower development in P. amurense. Two Cyc2CL transcripts (Cyc2CL-1 and Cyc2CL-2) regulate petal development and stamen abortion and are important for the ray floret (without stamens) development in the chrysanthemum [211]. RL-like gene, also named RAD, was one of the shared genes in both pathways leading to reversion from male to hermaphrodite flowers in hexaploid D. kaki and ectopic overexpression of DkRAD in model plants resulting in hypergrowth of the gynoecium [232]. PROLIFERA is necessary for megagametophyte and embryo development [233]. FBP2, a MADS-box gene involved in flower development, represents the same E function as SEP3 in Arabidopsis [234]. It has been demonstrated that YABBY genes participate in ovule development. Among YABBY genes, YABBY2 plays a critical role in ovule development and is expressed in different ovule development stages, while YABBY1 has the highest expression in the Megaspore Mother Cell stage [101]. It was reported that FIL, AS1, PHB, and YABBY1 interact with ANT and SEUSS to regulate ovule development and promote organ polarity. For example, ANT combined with the YABBY gene FIL to promote organ polarity by up-regulating the expression of the adaxial-specifying HD-ZIP gene PHABULOSA. The SEUSS/LUG coregulator complex physically interacts with YABBY1 and/or ANT to regulate adaxial identity genes PHB and REV to promote ovule and carpel growth [223][224][225][226]235]. Therefore, the genes above might form a flower development network, with ANT and SEUSS as initiation factors in P. amurense. Herein, a simple model of female flower development in P. amurense is proposed as a reference for future research (Figure 7). Additionally, FMO3, ATPD, WRKY21, and exostosin are located in the sex-determined region of plants with unknown specific functions [148,176,214]. Their specific functions need further study. Male-related genes TOZ19 and SERK1 were up-regulated in female and induced female flowers compared to male flowers [45,207]. The genes have not been well studied and whether they have the same role in P. amurense needs further verification. Furthermore, some male-related genes, including GPAT3, ZAT3, and FD [30, 83,89], were down-regulated, which might promote female development by suppressing male development.
Biochemical pathway-related genes and hormone-related genes might participate in each stage of female flower development in P. amurense. The two types of genes might contribute to extensive physiological activities. Their specific roles in P. amurense have not yet been well illustrated in this study.
Fan et al. discovered monoecious P. amurense for the first time and preliminarily explored the sexual formation and differentiation mechanism [30]. The discovery of monoecious P. amurense led to the sexual differentiation mechanism of P. amurense at the epigenetic level. This suggested that we should carefully observe the individual characteristics of plants, from which we can obtain inspiration. In this study, a precise method of demethylation was applied to further demonstrate that epigenetics was indeed involved in the sex regulation of P. amurense, and DNA methylation was one of its regulation modes.
Sex transition from male to female flowers with complete structure and function through demethylation suggested that demethylation might influence key genes for the sex differentiation in P. amurense. Among the genes analyzed, PAL3 and HUA1 were reported to regulate plant sex by DNA methylation [8,221]. Indeed, the performance of PAL3 and HUA1 in P. amurense after demethylation was consistent with the function reported. Therefore, these two genes deserve attention in subsequent studies.

Conclusions
Male flowers of P. amurense could be induced to convert into female flowers with complete structure and function with 5-azaC treatment, demonstrating that epigenetics is involved in the sex regulation of P. amurense. Potential sex-related genes were acquired from the transcriptome of P. amurense by referring to sex-related genes reported in the literature. Sex-related DEGs between female flowers (including F and F i ) and male flowers were identified by comparing female and male flowers. Further, 188 candidate genes of female flower development were screened out, including sex-regulatory genes, genes related to the formation and development of sexual organs, genes related to biochemical pathways, and hormone-related genes. Combined with gene functions and qRT-PCR analysis, RPP0W, PAL3, MCM2, MCM6, SUP, ANT, AIL6, AGL11, SEUSS, SRS5, and ESR2 were preliminarily considered to be the key genes for female flower development.